home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Technotools
/
Technotools (Chestnut CD-ROM)(1993).ISO
/
lang_c
/
cug167
/
spline.c
< prev
next >
Wrap
Text File
|
1985-08-21
|
15KB
|
658 lines
/*
HEADER: CUG
TITLE: SPLINE.C - Interpolate Smooth Curve
VERSION: 1.00
DATE: 11/01/85
DESCRIPTION: "SPLINE takes pairs of numbers from the standard
input as abscissae and ordinates of a function.
It produces a similar set, which is approximately
equally spaced and includes the input set, on the
standard output. The cubic spline output has two
continuous derivatives and sufficiently many
points to look smooth when plotted."
KEYWORDS: spline, graphics, plot, filter, UNIX
SYSTEM: ANY
FILENAME: SPLINE.DOC
WARNINGS: NONE
CRC: xxxx
SEE-ALSO: SPLINE.C
AUTHORS: Ian Ashdown - byHeart Software
COMPILERS: ANY
REFERENCES: AUTHORS: Bell Laboratories
TITLE: UNIX Programmer's Manual, Volume 1
p. 145; Holt, Rinehart and Winston
AUTHORS: R.W. Hamming
TITLE: Numerical Methods for Scientists and
Engineers, 2nd Edition
pp. 349 - 355; McGraw-Hill (1973)
AUTHORS: Forsythe, Malcolm & Moler
TITLE: Computer Methods for Mathematical
Computations
pp. 70 - 83; Prentice-Hall
ENDREF
*/
/*-------------------------------------------------------------*/
/* SPLINE.C - Interpolate Smooth Curve
*
* Version 1.00 November 1st, 1985
*
* Copyright 1985: Ian Ashdown
* byHeart Software
* 1089 West 21st Street
* North Vancouver, B.C. V7P 2C6
* Canada
*
* This program may be copied for personal, non-commercial use
* only, provided that the above copyright notice is included in
* all copies of the source code. Copying for any other use
* without previously obtaining the written permission of the
* author is prohibited.
*
* Synopsis: SPLINE [option] ...
*
* Description: SPLINE takes pairs of numbers from the standard
* input as abscissae and ordinates of a function.
* It produces a similar set, which is approximately
* equally spaced and includes the input set, on the
* standard output. The cubic spline output (R.W.
* Hamming, "Numerical Methods for Scientists and
* Engineers", 2nd ed. 349ff) has two continuous
* derivatives and sufficiently many points to look
* smooth when plotted.
*
* The following options are recognized, each as a
* separate argument:
*
* -a Supply abscissae automatically (they are
* missing from the input); spacing is given by
* the next argument, or is assumed to be 1 if
* next argument is not a number.
*
* -k The constant "k" is used in the boundary
* value computation
*
* y " = ky " , y " = ky "
* 0 1 n n-1
*
* is set by the next argument. By default,
* k = 0.
*
* -n Space output points so that approximately n
* intervals occur between the lower and upper
* "x" limits (default n = 100.)
*
* -p Make output periodic, i.e. match derivatives
* at ends. First and last input values must
* agree.
*
* -x Next 1 (or 2) arguments are lower (and upper)
* "x" limits. Normally these limits are
* calculated from the data. Automatic abscissae
* start at lower limit (default 0).
*
* Diagnostics: When data is not strictly monotone in "x", SPLINE
* reproduces the input without interpolating extra
* points.
*
* Bugs: A limit of 1000 input points is silently
* enforced.
*
* (Note: the above description is a slightly reworded version of
* that appearing in the "UNIX Programmer's Manual",
* copyright 1979, 1983 Bell Laboratories.)
*/
/*** Definitions ***/
#define FALSE 0
#define TRUE 1
#define MAX_SIZE 1000 /* Input point array limit */
#define ILL_ARG 0 /* Error codes */
#define ILL_CMB 1
#define ILL_KVL 2
#define ILL_OPT 3
#define INS_INP 4
#define MIS_KVL 5
#define MIS_XVL 6
#define MIS_YVL 7
#define NMT_ORD 8
#define SQUARE(a) a*a
#define CUBE(a) a*a*a
/*** Typedefs ***/
typedef int BOOL; /* Boolean flag */
/*** Include Files ***/
#include <stdio.h>
#include <ctype.h>
#include <math.h>
/*** Main Body of Program ***/
int main(argc,argv)
int argc;
char **argv;
{
int n = 0,
i,
j,
n_val = 0;
float x[MAX_SIZE],
y[MAX_SIZE];
double a_val = 1.0,
k_val = 0.0,
x1_val = 0.0,
x2_val = 0.0,
x_intvl,
ix,
iy,
d2y[MAX_SIZE],
h,
atof(),
fabs(),
spl_int();
char buffer[257],
*temp,
*gets();
BOOL aflag = FALSE, /* Command-line option flags */
kflag = FALSE,
pflag = FALSE,
x1flag = FALSE,
x2flag = FALSE,
is_float();
void spl_coeff(),
pspl_coeff(),
error();
/* Parse the command line for user-selected options */
while(--argc)
{
temp = *++argv;
if(*temp != '-') /* Check for legal option flag */
error(ILL_OPT,*argv);
else
switch(toupper(*++temp))
{
case 'A': /* "-a" option */
aflag = TRUE;
if(argc > 1 && is_float(*(argv+1)))
{
argc--;
argv++;
if((a_val = atof(*argv)) < TINY)
error(ILL_ARG,*argv);
}
break;
case 'K': /* "-k" option */
if(pflag == TRUE)
error(ILL_CMB,NULL);
kflag = TRUE;
if(argc > 1 && is_float(*(argv+1)))
{
argc--;
argv++;
k_val = atof(*argv);
if(k_val - 4.0 < 0.000001)
error(ILL_KVL,*argv);
break;
}
else
error(MIS_KVL,NULL);
case 'P': /* "-p" option */
if(kflag == TRUE)
error(ILL_CMB,NULL);
pflag = TRUE;
break;
case 'X': /* "-x" option */
x1flag = TRUE;
if(argc > 1 && is_float(*(argv+1)))
{
argc--;
argv++;
x1_val = atof(*argv);
}
else
error(MIS_XVL,NULL);
if(argc > 1 && is_float(*(argv+1)))
{
x2flag = TRUE;
argc--;
argv++;
x2_val = atof(*argv);
}
break;
default: /* "-n" option */
while(isdigit(*temp))
n_val = n_val * 10 + (*temp++ - '0');
if(*temp != '\0')
error(ILL_OPT,*argv);
}
}
if(n_val == 0) /* Set "n_val" if not given */
n_val = 100;
/* Get the input data */
while(1) /* ... while there is more input data */
{
if(aflag == TRUE) /* Automatic abscissae were called for */
{
if(n == 0)
x[0] = x1_val;
else
x[n] = x[n-1] + a_val;
}
else /* Abscissae supplied with input data */
{
if(gets(buffer))
x[n] = atof(buffer);
else
break;
}
if(gets(buffer)) /* Read in the corresponding ordinate */
y[n] = atof(buffer);
else
{
if(aflag == TRUE)
break;
else
error(MIS_YVL,NULL);
}
if(++n == MAX_SIZE) /* Maximum amount of input data? */
break;
}
if(n <= 2) /* Check for insufficient input data */
error(INS_INP,NULL);
/* Check for non-monotonic abscissae. Output input data set
* without interpolation if true.
*/
h = x[1] - x[0];
for(i = 1; i < n-1; i++)
if(fabs(x[i+1] - x[i] - h) > 0.0)
{
for(i = 0; i < n; i++)
printf("%g\n%g\n",x[i],y[i]);
exit();
}
/* Calculate abscissa interval. Use "-x" option values if
* they were given.
*/
if(x1flag == FALSE)
x1_val = x[0];
if(x2flag == FALSE)
x2_val = x[n-1];
x_intvl = (x2_val - x1_val)/(n_val - 1);
/* Find the coefficients */
if(pflag == FALSE)
spl_coeff(y,d2y,h,n,k_val);
else
pspl_coeff(y,d2y,h,n);
/* Interpolate and output results */
ix = x1_val;
i = 1;
for(j = 0; j < n_val; j++)
{
if(ix >= x[i] && i < n - 1)
i++;
iy = spl_int(ix,x,y,d2y,h,i);
printf("%g\n%g\n",ix,iy);
ix += x_intvl;
}
}
/*** Functions ***/
/* SPL_COEFF() - Calculate spline coefficients and return in
* vector "d2y[]". Matrix to be solved has the
* typical form:
*
* +- -+ +- -+ +- -+
* | 4-k 1 0 0 | | y1" | = m * | y2 - 2*y1 + y0 |
* | 1 4 1 0 | | y2" | | y3 - 2*y2 + y1 |
* | 0 1 4 0 | | y3" | | y4 - 2*y3 + y2 |
* | 0 0 1 4-k | | y4" | | y5 - 2*y4 + y3 |
* +- -+ +- -+ +- -+
*
* where k = k_val, m = 6.0/(h*h) and yn" is the
* second derivative of the interpolated function
* at the "nth" abscissa, y0" = k * y1" and
* y5" = k * y4".
*/
void spl_coeff(y,d2y,h,n,k_val)
float y[];
double d2y[],
h,
k_val;
int n;
{
double m,
a[MAX_SIZE-1];
int i;
/* Set up the (symmetric tridiagonal) matrix, where the only
* elements of interest are those on the diagonal. These are
* stored in array "a[]". Array "d2y[]" initially holds the
* constants vector, then is overlaid with the calculated
* variables vector.
*/
m = 6.0/(h*h);
for(i = 1; i < n-1; i++)
d2y[i] = m * (y[i+1] - 2.0 * y[i] + y[i-1]);
a[1] = 4.0 - k_val;
/* Reduce the matrix to upper triangular form */
for(i = 2; i < n-2; i++)
{
a[i] = 4.0 - 1.0/a[i-1];
d2y[i] -= d2y[i-1]/a[i-1];
}
a[n-2] = 4.0 - k_val - 1.0/a[n-3];
d2y[n-2] -= d2y[n-3]/a[n-3];
d2y[n-2] /= a[n-2];
/* Solve using back substitution */
for(i = n-3; i > 0; i--)
d2y[i] = (d2y[i] - d2y[i+1])/a[i];
/* Solve for end conditions */
d2y[n-1] = d2y[n-2] * k_val;
d2y[0] = d2y[1] * k_val;
}
/* PSPL_COEFF() - Calculate periodic spline coefficients and
* return in vector "d2y[]". Matrix to be solved
* has the typical form:
*
* +- -+ +- -+ +- -+
* | 4 1 0 0 1 | | y1" | = 6 * | y2 - 2*y1 + y0 |
* | 1 4 1 0 0 | | y2" | | y3 - 2*y2 + y1 |
* | 0 1 4 1 0 | | y3" | | y4 - 2*y3 + y2 |
* | 0 0 1 4 1 | | y4" | | y5 - 2*y4 + y3 |
* | 1 0 0 1 4 | | y5" | | y1 - y0 - y5 + y4 |
* +- -+ +- -+ +- -+
*
* where m = 6.0/(h*h) and yn" is the second
* derivative of the interpolated function at the
* "nth" abscissa and y0" = y5".
*/
void pspl_coeff(y,d2y,h,n)
float y[];
double d2y[],
h,
fabs();
int n;
{
double c,
m;
static double a[MAX_SIZE-1],
b[MAX_SIZE];
int i;
/* Check for matching end ordinates */
if(fabs(y[n-1] - y[0]) > 0.0)
error(NMT_ORD,NULL);
/* Set up the matrix, where array "a[]" holds the diagonal
* elements, "b[]" holds those elements of column "n-1", and
* "c" holds the current element of interest of row "n-1".
* Array "d2y[]" initially holds the constants vector, then is
* overlaid with the calculated variables vector.
*/
d2y[0] = 0.0;
m = 6.0/(h*h);
for(i = 1; i < n-1; i++)
d2y[i] = m * (y[i+1] - 2.0 * y[i] + y[i-1]);
d2y[n-1] = m * (y[1] - y[0] - y[n-1] + y[n-2]);
a[1] = 4.0;
b[1] = 1.0;
b[n-2] = 1.0;
b[n-1] = 4.0;
c = 1.0;
/* Reduce the matrix to upper triangular form */
for(i = 2; i < n-1; i++)
{
m = 1/a[i-1];
a[i] = 4.0 - m;
b[i] -= b[i-1] * m;
d2y[i] -= d2y[i-1] * m;
b[n-1] -= c * m * b[i-1];
d2y[n-1] -= c * m * d2y[i-1];
c = -c * m;
}
c += 1.0;
b[n-1] -= c * b[n-2]/a[n-2];
d2y[n-1] -= c * d2y[n-2]/a[n-2];
/* Solve using back substitution */
d2y[0] = d2y[n-1] /= b[n-1];
d2y[n-2] = (d2y[n-2] - b[n-2] * d2y[n-1])/a[n-2];
for(i = n-3; i > 0; i--)
d2y[i] = (d2y[i] - b[i] * d2y[n-1] - d2y[i+1])/a[i];
}
/* SPL_INT - Interpolate points using spline function */
double spl_int(ix,x,y,d2y,h,i)
float x[],
y[];
double ix,
d2y[],
h;
int i;
{
double iy,
t1,
t2;
t1 = (ix - x[i-1])/h;
t2 = (x[i] - ix)/h;
iy = y[i-1] * t2 + y[i] * t1 - SQUARE(h) * (d2y[i-1] * (t2 -
CUBE(t2)) + d2y[i] * (t1 - CUBE(t1)))/6.0;
return iy;
}
/* IS_FLOAT() - Check that character string is in correct floating
* point format. Return TRUE if correct, FALSE
* otherwise. The algorithm used is a deterministic
* finite state machine. Using the regular
* expression terminology of Unix's "lex", the
* character string must be of the form:
*
* -?d*.?d*(e|E(\+|-)?d+)?
*
* where d = 0|1|2|3|4|5|6|7|8|9
*/
BOOL is_float(str)
char *str;
{
int c, /* Next FSM input character */
state = 0; /* Current FSM state */
while(c = *str++)
{
switch(state)
{
case 0: /* FSM State 0 */
switch(c)
{
case '-':
state = 1;
break;
case '.':
state = 3;
break;
default:
if(isdigit(c))
state = 2;
else
return FALSE;
break;
}
break;
case 1: /* FSM State 1 */
switch(c)
{
case '.':
state = 2;
break;
default:
if(isdigit(c))
state = 2;
else
return FALSE;
break;
}
break;
case 2: /* FSM State 2 */
switch(c)
{
case '.':
state = 4;
break;
case 'e':
case 'E':
state = 5;
break;
default:
if(isdigit(c))
state = 2;
else
return FALSE;
break;
}
break;
case 3: /* FSM State 3 */
if(isdigit(c))
state = 4;
else
return FALSE;
break;
case 4: /* FSM State 4 */
switch(c)
{
case 'e':
case 'E':
state = 5;
break;
default:
if(isdigit(c))
state = 4;
else
return FALSE;
break;
}
break;
case 5: /* FSM State 5 */
switch(c)
{
case '+':
case '-':
state = 6;
break;
default:
if(isdigit(c))
state = 7;
else
return FALSE;
break;
}
break;
case 6: /* FSM State 6 */
if(isdigit(c))
state = 7;
else
return FALSE;
break;
case 7: /* FSM State 7 */
if(isdigit(c))
state = 7;
else
return FALSE;
break;
}
}
return TRUE;
}
/* ERROR() - Error reporting. Returns an exit status of 2 to the
* parent process.
*/
void error(n,str)
int n;
char *str;
{
fprintf(stderr,"\007\n*** ERROR - ");
switch(n)
{
case ILL_ARG:
fprintf(stderr,"Argument must be greater than zero: %s",
str);
break;
case ILL_CMB:
fprintf(stderr,"Cannot use -k option with -p option");
break;
case ILL_KVL:
fprintf(stderr,"Value for -k option cannot equal 4.0: %s",
str);
break;
case ILL_OPT:
fprintf(stderr,"Illegal command line option: %s",str);
break;
case INS_INP:
fprintf(stderr,"Insufficient input data");
break;
case MIS_KVL:
fprintf(stderr,"Missing value for -k option");
break;
case MIS_XVL:
fprintf(stderr,"Missing value for -x option");
break;
case MIS_YVL:
fprintf(stderr,"Missing ordinate value");
break;
case NMT_ORD:
fprintf(stderr,"End ordinates do not match");
break;
default:
break;
}
fprintf(stderr," ***\n\nUsage: spline [-aknpx]\n");
exit(2);
}
/*** End of SPLINE.C ***/
fprintf(stderr,"\007\n*** ERROR - ");
switch(n)
{
case ILL_ARG:
fprintf(stderr,